Models for diffusion and island growth in metal monolayers 
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A model that describes self diffusion, island nucleation and film growth on FCC(OOf) metal 
substrates is presented. The parameters of the model are optimized to describe Cu diffusion on 
Cu(OOf), by comparing activation energy barriers to a full set of barriers obtained from semi- 
empirical potentials via the embedded atom method. It is found that this model (model I), with 
only three parameters, provides a very good description of the full landscape of hopping energy 
barriers. These energy barriers are grouped in four main peaks. A reduced model (model II) with 
only two parameters, is also presented, in which each peak is collapsed into a single energy value. 
From the results of our simulations, we find that this model still maintains the essential features of 
diffusion and growth on this model surface. We find that hopping rates along island edges are much 
higher than for isolated atoms (giving rise to compact island shapes) and that vacancy mobility 
is higher than adatom mobility. We observe substantial dimer mobility (comparable to the single 
atom mobility) as well as some mobility of trimers. Mobility of small islands affects the scaling 
of island density N vs. deposition rate F, N ~ F 1 , as well as the island size distribution. In the 
asymptotic limit of slow deposition, scaling arguments and rate equations show that 7 = i* /(2i* + 1) 
where i* is the size of the largest mobile island. Our Monte Carlo results, obtained for a range of 
experimentally relevant conditions, show 7 = 0.32 ± 0.01 for the EAM, 0.33 ± 0.01 for model I and 
0.31 ± 0.01 for model II barriers. These results are lower than the anticipated 7 > 0.4 due to dimer 
(and trimer) mobility. 
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I. INTRODUCTION 

Recent experiments on thin film growth on well char- 
acterized substrates using molecular beam epitaxy have 
provided detailed information about growth kinetics and 
morphology. In particular, growth in the submonolayer 
regime has been studied extensively using both scanning 
tunneling microscopy (STM) and diffraction meth- 

ods such as helium atom beam scattering jll]-|l3) and low 
energy electron diffraction 14 Tq|. It was observed that 
for a variety of systems and a broad temperature range, 
island nucleation is the dominant mechanism for crystal 
growth. A variety of island morphologies has been found. 
Fractal- like islands, resembling diffusion limited aggrega- 
tion (DLA) clusters, were observed in Au on Ru(0001) 
|,§, Cu on Ru(0001) § and Pt on Pt(lll) §§. On 
the other hand, compact islands were observed in Ni on 
Ni(OOl) system ||. Scaling properties of the island den- 
sity as a function of deposition rate and coverage, as well 
as the island size distribution have be en stu died exper- 
imentally by several groups [0,|],[lO [ 
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microscopy (FIM) experiments that provide direct access 
to diffusion processes at the atomic scale, have also been 
performed fl7|-|l9||. This technique was used to identify 
the diffusion modes of adatoms |l^] as well as small is- 
lands |ll|[ll| on FCC(OOl) metal surfaces and to measure 
their diffusion coefficients. 

Theoretical studies aimed at providing a better un- 
derstanding of the relation between key processes at the 
atomic scale and the resulting morphologies have been 
done using Monte Carlo (MC) simulations [po|^5|. In 
simulations of island growth during deposition, atoms are 
deposited randomly on the substrate at rate F [given in 
monolayer (ML) per second] and then hop according to a 
microscopic model. The hopping rate h (in units of hops 
per second) of a given atom to each unoccupied nearest 
neighbor (NN) site is given by 



h = v ■ cxp(— Es/ksT) 



(1) 



where v = 10 12 s is the commonly used attempt rate, 
Eb is the activation energy barrier, ks is the Boltzmann 
constant and T is the temperature. 

The activation energy barrier Eb depends on the lo- 
cal environment of the hopping atom, namely the con- 
figuration of occupied and unoccupied adjacent sites. 
Two approaches have been taken in the construction of 
the energy barriers for hopping. One approach was to 
construct simple models that include the desired fea- 
tures, such as stability and mobility of small islands, 
and that take into account properties such as bond en- 
ergies @-||,||,||,|||l|,|| |§. The advantage of this 
approach is that the models are well defined and use only 
a few parameters. These models are useful for studies of 
scaling and morphology but cannot provide a quantita- 
tive description of diffusion on a particular substrate. A 
second approach is based on the use of an approximate 
many-body energy functional to calculate the hopping 



energy barriers for a complete set of relevant configura- 
tions |24|,^],!2| . This approach provides a good descrip- 
tion of diffusion processes on the given substrate but only 
limited understanding due to the large number of param- 
eters. 

In this work we propose a procedure that combines 
the advantages of both approaches. Using sensible as- 
sumptions about the bond energies and diffusion paths 
we obtain a simple formula for the activation energy bar- 
riers. We then optimize the parameters of this formula by 
using energy barriers obtained from the embedded atom 
method for Cu diffusion on Cu(OOl). This procedure pro- 
vides two simple models which combine the best features 
of both approaches. Model I, which has three parameters, 
provides good quantitative description of the landscape 
of hopping energy barriers. Model II, which has only two 
parameters and is as simple as other minimal models, still 
incorporates the essential physics of diffusion of adatoms 
on FCC(OOl) metal surfaces. 

We find that on the Cu(OOl) substrate, edge mobility 
(i.e., the mobility of an adatom at the perimeter of an 
island) is much higher than the single adatom mobility, 
giving rise to compact island shapes. Dimer mobility is 
comparable to the single adatom mobility while trimers 
are also mobile [ p7||3^j37[ | . Mobility of small islands has 
a significant effect on the asymptotic scaling properties 
and is thus particularly important. We also find that 
vacancies have higher mobility than adatoms. 

The paper is organized as follows. In Section II we in- 
troduce the models. Simulation and results are presented 
in Section III with emphasis on scaling and morphology. 
The results and their implications are discussed in Sec- 
tion IV and a summary is given in Section V. 



II. THE MODELS 
A. The EAM Barriers 

In this work we use a set of energy barriers for Cu on 
Cu(OOl) |38|l obtained using the embedded atom method 
(EAM) |39| . This method uses semiempirical potentials 
and provides a good description of self diffusion of Cu 
on Cu(OOl) and similar surfaces pq |. Specifically, we use 
the EAM functions of Cu developed by Adams, Foiles, 
and Wolfer |^] which are fitted to a similar data base as 
the one employed by Foiles, Baskes, and Daw EJ . 

The hopping energy barriers are calculated for all local 
environments as shown in Fig. 1, where seven adjacent 
sites, i = 0, . . . , 6 are taken into account |l2|]. Each one 
of these sites can be either occupied (Si = 1) or vacant 
(Si = 0), giving rise to 2 7 = 128 barriers. A binary 
representation is used to assign indices to these barriers. 
For each configuration (So, . . . , Sq) the barrier is given 



by Eg, where 
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n = £ Si ■ T (2) 

takes the values n — 0, . . . , 127. The full set of hop- 
ping energy barriers (given in eV) is presented in Table 
I. To show these values in a compact form, each barrier 
in Table I corresponds to a configuration in which the 
occupied sites are the union of the occupied sites in the 
picture on top of the given column and on the left hand 
side of the given row. The column in Table I in which 
a given configuration appears is determined by the occu- 
pancy of sites i = 2, 3, 6 while the row is determined by 
sites i = 0, 1, 4, 5. One can define 

i=2,3,6 1=0,1,4,5 

such that for each configuration n = ri\ +ri2- To demon- 
strate the use of Table I, we will check the barrier for the 
configuration in which sites 0, 3 and 4 are occupied and 
all other sites adjacent to the hopping atom are vacant. 
For this configuration, according to Eq. (||), n\ = 8 and 
ii2 = 17 [n = 25). The barrier, which is found in the 
column with the index 8 and the row with index 17, is 
Ef = 0.89 eV. 

In Table I we use the symmetries of the configurations 
in the 3x3 cell to reduce the number of entries. There is a 
mirror symmetry plane perpendicular to the surface and 
containing the arrow of the hopping atom. Consequently, 
the columns of n\ = 4 and 12, in which site i = 2 is 
occupied, stand also for the symmetric configurations in 
which i — 6 is occupied. In the other four columns, there 
are some configurations that, due to symmetry, appear 
twice. In such cases, the barrier for the configuration 
with larger n appears in italics. 

Here we consider only hopping moves in which a single 
atom hops each time. However, it turns out that in some 
cases the molecular statics calculations, used to obtain 
the barriers, give rise to concerted moves. In such moves 
the atom at site i — 3 follows the hopping atom and 
takes the place vacated by the hopping atom. This fact 
significantly reduces the barrier. We have found that for 
all configurations in which concerted moves appear, these 
concerted moves can be suppressed by adding one more 
row of atoms on the left hand sides of sites i = 0, 3 and 
4. In Table I, the energy values for those configuration in 
which a concerted move was found, are shown in paren- 
thesis. The barrier obtained when the concerted move 
was suppressed is shown to the left of the parenthesis. 

Concerted moves may affect film growth by modifying 
either the stability and mobility of small islands or the 
morphology of large islands. The effect of the shear move 
(ni = 8; ri2 = 3), in which concerted motion reduces 
the barrier from E^ 1 = 0.78 eV to 0.65 eV, was studied 
in Ref. |43] ]. It was claimed that for experiments with 
very low deposition rates this move effectively raises i* 
from 3 to 8, where i* is the size of the critical island 
nucleus. Simple calculations of the hopping rates, using 



Eq. (|l|) indicate that, since the reduced barrier is still 
rather large, the effect of this concerted move is negligible 
for the temperatures and deposition rates studied in this 
paper . 

We find that there are other concerted moves which, 
unlike the shear move discussed above, cause a dramatic 
reduction in the barrier for the given move (Table I). 
If the atoms inside the 3x3 square represent an iso- 
lated island and sites around are vacant, it is easy to see 
that moves such as (ni,n 2 ) = (34, 12), and (35,76) sim- 
ply lead to a more stable configuration of the pentamer 
and heptamer islands, respectively. Once this more sta- 
ble configuration is obtained, the concerted move cannot 
take place again. Since these more stable configurations 
can be obtained through other fast moves, the concerted 
move does not have a significant effect on the overall mor- 
phology of the film. Other moves such as (34, 8), (50, 12), 
and (51, 76) can occur only as long as the islands of 4, 6 
and 8 atoms respectively, have not reached their most sta- 
ble configuration. If the above configurations, in which 
concerted moves are possible, occur in a denser environ- 
ment, and in particular if the sites on the left hand side 
of sites i = 0, 3 and 4 are occupied, the concerted moves 
are suppressed. Thus, we conclude that in both cases 
the concerted moves do not have a significant effect in 
the simulations presented here and they will be ignored. 

To gain a better understanding of the barrier en- 
ergy landscape we discuss the barrier height distribution 
(without concerted moves) in Fig. 2(a). We observe that 
this distribution exhibits four peaks. This feature is in 
agreement with Ref. [Q where a different method flsj ] 
was used to calculate the barriers. Each peak, corre- 
sponds to a single or a double column in Table I (there is 
a little ambiguity in the region between peaks III and IV 
due to a slight overlap). In general, peak I includes very 
fast moves toward island edges, peak II includes moves 
along the edge, peak III includes, most notably, the single 
atom move , while peak IV includes detachment moves. 

B. Construction of Models 

When an atom on the surface hops into a vacant near- 
est neighbor site it has to cross the energy barrier be- 
tween the initial and final sites. We have used molecular 
statics in conjunction with EAM functions to find the 
diffusion path of an adatom. Not surprisingly, we have 
found that the top of the energy barrier to go from a four- 
fold coordinate site to an adjacent site is at the bridge 
site. By slowly moving the adatom within the molecular 
statics calculation, it turns out that the hopping energy 
barrier is simply the difference between the energy at the 
bridge site and in the initial site. The occupancy of the 
seven adjacent sites (Fig. 1) affects both energies. We 
will now express the energy of the hopping atom in its 
initial site as: 
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E in — E in — AE in ■ (Si + S3 + S5) 

-AE™ in ■ (So + S 2 + S± + S 6 ) (4) 

where Si = 1 if site i is occupied and otherwise. The 
energy of an isolated atom is E® n while the reduction in 
its energy due the presence of an atom in a nearest (next 
nearest) neighbor site is given by AE in (AEV^ n ). Here 
we assume that the contributions of nearest neighbor and 
next nearest neighbor (NNN) atoms to the energy are 
additive. The energy of the hopping atom when it is on 
the bridge site is given by: 

Etop = E® op — AE top ■ (Si + S2 + S5 + 5*6) (5) 

where E^ op is the energy of an isolated atom on top of 
a bridge site, while AE top is the reduction in the energy 
due to the presence of an atom in one of the four sites 
adjacent to the bridge site. We do not include here NNN 
type contributions since their effect is small. Therefore, 
for a given configuration the barrier Eb — E top ~ Ei n for 
an atom to hop into an adjacent vacant site is given by: 

E 7 g = E B — AE to p ■ (S2 + So) + AE in ■ S3 
— (AE top — AE in ) ■ (Si + S5) 
+AEZ n ■ (So + S 2 + S 4 + S 6 ) (6) 

where E B = E^ op — Ef n and n is given by Eq. (||). To 
examine the formula above, we found the parameters that 
best describe the 128 EAM barriers by minimizing the 
sum of squares 

127 

R = Y J [E n B {EAM)-El(Eq. 6)] 2 . (7) 

n=0 

The values obtained for the parameters in Eq. (||) 
are E% = 0.494, AE in = 0.265, AE top = 0.268 and 
AE™ n = 0.024 eV. Thus, we find that to within about 
0.003 eV, AE in = AE top . Replacing both AE m and 
AE top by AE — 0.265 eV we obtain a three parameter 
model (model I) in which the energy barriers are: 

Eb — E B + AE ■ {S 3 -S 2 - S 6 ) 

+AE™ in ■ (So + S 2 + S 4 + S 6 ) (8) 

where n is given by Eq. (^|). The distribution of energy 
barriers obtained from Eq. (||) is shown in Fig. 2(b). 
A very good agreement in the location of the four peaks 
with the EAM energy barriers is found, but EAM peaks 
are significantly broader. This can be due to the fact 
that during the hopping move adjacent atoms can relax 
within their potential well. Model I accounts for these ef- 
fects only on average and therefore gives rise to narrower 
peaks. 

One can further simplify the model by ignoring the 
NNN interactions which are relatively small, namely 
choosing AEl" ln = 0. This model (model II) has only 
two parameters and is described by: 



Eb = E B + AE ■ (S 3 - S 2 - S 6 ) (9) 

where n is given by Eq. (Q). Repeating the optimization 
procedure of Eq. (f7j) for the barriers given by Eq. (|]), 
we obtain E B = 0.526 and AE = 0.255. In this model, 
all the energy barriers in each peak collapse into a single 
value. In spite of its simplicity, this model provides a 
good description of self diffusion of Cu on the Cu(001) 
substrate. 



C. Diffusion Coefficients 

To find out which island sizes are mobile and to obtain 
the diffusion coefficients of mobile islands on Cu(001), 
we have done simulations of single cluster diffusion. To 
obtain the statistics required for a precise determina- 
tion of the diffusion coefficients we performed 1000 runs 
for diffusion of monomers, dimers, trimers, pentamers 
(,s = 5) and hcptamers (s = 7). Each run was carried 
out for a time equal to 1.0 seconds, which is more than 
100 times larger than the time scale for hopping of an 
isolated adatom at the given temperature (T — 250K). 
The diffusion coefficients were obtained from the relation 
(r 2 ) = 4h s t, where s — 1, 2, 3, 5, 7 is the number of atoms 
in the cluster, r is the distance between the initial posi- 
tion of the center of mass of the cluster and its position 
after time t, h s is the diffusion coefficient for a cluster of 
size s and (...) represents an average over the 1000 runs. 
The diffusion coefficients for monomers, dimers, trimers, 
pentamers and heptamers are shown in Table II. It is also 
found that under the conditions studied here the diffu- 
sion coefficients for islands of size s — 4, 6 and s > 8 are 
negligible. 

III. SIMULATIONS AND RESULTS 
A. Monte Carlo Simulations 

To assess the validity of our models, we have performed 
MC simulations of island growth for a range of deposition 
rates using the EAM as well as model I and model II 
barrier s |j46[| . We used a continuous time MC technique 
|g^-|^,|7[in which moves are selected randomly from 
the list of all possible moves at the given time with the 
appropriate weights. The time is advanced after each 
move according to the inverse of the sum of all rates. 
The existence of four peaks in the spectrum of energy 
barriers indicates that there are four typical time scales 
of hopping. The two lowest peaks include very fast moves 
towards and along island edges and motion of vacancies. 
The single atom move is in the third peak while moves in 
which atoms detach from islands are in the highest peak. 

From statistics collected during the simulations we find 
that most of the computer time is consumed by moves 
n =6 (96) and 7 (112). The first move occurs in trimers 
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while both moves occur for atoms hopping along straight 
island edges. The reason that these moves consume so 
much time is that the reverse move typically has the same 
low barrier so the atom can continuously hop back and 
forth for a long time. To make the simulation feasible, 
we had to somewhat suppress these moves by artificially 
raising the barriers in the first and second peaks. This 
was done for all barriers lower than 0.4 eV according 
to: Eg — > Eg + a(0A - E%) where a = 0.7. Since 
the moves associated with these barriers are still orders 
of magnitude faster than for the two higher peaks we 
expect that this modification will have only small effect 
on the island morphology. We tried various values of a 
and found that up to a » 0.8 this is indeed the case. 



barriers for three deposition rates, F = 0.1, 0.01 and 
0.001 ML/s are presented in Fig. 5(a)-(c), respectively. 
The simulations were done on a 250 x 250 lattice at 
T = 250K and the snapshots were taken at 8 = 0.1 
ML. It is observed that decreasing the deposition rate 
gives rise to fewer and larger islands. In Fig. 6 we 
present the island density vs. coverage for the EAM 
barriers, at deposition rates (from top to bottom) F = 
0.1,0.04,0.01,0.004,0.001 and 0.0004 ML/s. It is ob- 
served that in addition to the decrease in island den- 
sity as the deposition rate decreases, the plateau becomes 
broader and flatter. 



C. Scaling Properties 



B. Island Growth and Morphology 

In Fig. 3 we show the island morphologies obtained 
in MC simulations using the energy barriers from EAM 
[Fig. 3(a)], model I [Fig. 3(b)] and model II [Fig. 3(c)]. 
These simulations were done on a 250 x 250 lattice at 
T = 250K and deposition rate F = 0.01 ML/s. The 
snapshots presented here are for a coverage 8 = 0.1 ML. 
These snapshots indicate that models I and II provide a 
good description of the model with the full EAM barri- 
ers as far as the island morphology is concerned. The 
islands are rather compact, dominated by overall square 
and rectangular shapes with a small number of kinks. 

The island density vs. coverage is presented in Fig. 4 
for the EAM (*), model I (o) and model II (+) energy 
barriers. In all cases the island density quickly increases 
at low coverage, then saturates and remains nearly con- 
stant thereafter. It starts to decrease at higher coverage, 
when coalescence sets in. For model I the island density is 
slightly higher than for the EAM barriers while for model 
II it is about 60% higher. These differences are not intrin- 
sic to the models. They can be traced to the fact that 
the single atom hopping energy barriers obtained from 
the optimization procedure that we use are E B = 0.494 
eV for model I and E B = 0.526 eV for model II; these 
values are higher than the value of E B = 0.485 eV from 
the EAM calculations. These discrepancies could be fixed 
by an overall rescaling of the barriers in each model so 
that the single atom barrier would be exactly equal to 
the EAM value. Since there is some arbitrariness in this 
procedure we chose not to apply it in the simulations 
presented here. The important result of this work is that 
our models reproduce the essential features of diffusion, 
island growth and scaling. The overall factor in island 
density to make the models agree with the EAM calcu- 
lations could be obtained by rescaling. While model I 
reproduces the spectrum of EAM energy barriers more 
accurately, the importance of model II is that it captures 
the essential features of adatom diffusion on the Cu(001) 
in spite of its simplicity . 

Snapshots of the island morphology using the EAM 



The island density vs. deposition rate is presented in 
Fig. 7 for the EAM (*), model I (o) and model II (+) 
energy barriers for T = 250K and 8 = 0.1 ML. This 
coverage is well within the quasi-steady state regime in 
which the island density is nearly constant (see Fig. 6). 
In this regime scaling arguments and rate equations in- 
dicate that in the asymptotic limit of low deposition rate 
the island density exhibits power law behavior of the form 



N 



F 



(10) 



where h\ is the adatom diffusion coefficient (the hop- 
ping rate for an isolated atom is Ah\ due to the four 
possible directions for hopping). The exponent 7 de- 
pends on microsc opic p roperties of the system such as 
stability Jsij- lssl , !^ , ^ , ^ and mobility (5^j3^] of small is- 
lands, anisotropic diffusion |IpC[] and magic islands (is- 
lands which are stable while larger islands are unstable) 
p9[ . It was found that for systems in which the smallest 
stable island is of size i* + 1 (where islands of size s < i* 
dissociate), 7 = + 2). In this case i* is the critical 



island size 51 1. 

Under the conditions studied here, detachment of 
atoms from islands is negligible, while mobility of small 
islands is a significant effect. Therefore, we define here 
the critical island size i* as the size for which all islands 
of size s < i* are mobile while islands of size s > i* + 1 
are immobile. In this case the asymptotic value of 7 in 
the limit where F/hi — > is given by [p5 35 : 



7 



2r f 



1 



(11) 



This result is exact if all mobile islands have the same 
diffusion coefficient, namely h s = h, s = 1, . . . , i*. Still, 
it is a good approximation if all the diffusion coefficients 
are of the same order of magnitude. Table II shows that 
the diffusion coefficients of pentamers and heptamers are 
practically negligible. For EAM the diffusion coefficient 
of the trimer is also rather small and we expect the sys- 
tem to be described by i* = 2. Models I and II exhibit 
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small but not negligible mobility of trimers and there- 
fore we expect them to be described by either i* = 2 or 
i*=3. 

In Fig. 7 we observe nearly two and a half decades of 
scaling behavior, in a range of deposition rates between 
F = 0.1 and 0.0004 ML/s. The best fit through the 
six data points provides 7 = 0.32 ± 0.01 for the EAM, 
0.33 ± 0.01 for model I and 0.31 ± 0.01 for model II barri- 
ers. These results are significantly lower than the asymp- 
totic values for i* = 2 (7 = 2/5) or i* = 3 (7 = 3/7). 
In Ref. |35| a similar deviation was found between simu- 
lation results and the asymptotic predictions of the rate 
equations. It was argued that this deviation is due to the 
slow convergence of 7 to its asymptotic value, when the 
deposition rate F is lowered. This conclusion is strength- 
ened by the fact that if only the four left-most data points 
are used in the linear fit, the resulting 7 equals 0.35±0.01 
for the EAM, 0.36 ± 0.01 for model I and 0.34 ± 0.01 
for model II. These results are somewhat closer to the 
asymptotic value. Moreover, within the error bars they 
are larger than the asymptotic value that arises from 
adatom mobility alone, namely 7 = 1/3. Therefore, the 
possibility that 7 converges to 1/3 should be discarded. 
Since the conditions used here are experimentally rele- 
vant, our results indicate that one should be careful in 
drawing conclusions about processes at the atomic scale 
from scaling results (see Section IV D). 

The scaling properties of the island size distribution 
have been studied both experimentally 14 1§ and theo- 
retically j|2|||||Q|3,|||. These studies indicated that 
the island size distribution depends on the stability and 
mobility of small islands and is modified in the case of 
magic island sizes |2^]. The scaled island size distribu- 
tions are presented for the three models in Fig. 8. For 
all three models the deposition rates are F = 10 _1 (*), 
10~ 2 (o) and 10~ 3 ML/s (+). These results are based 
on statistics collected from 20 runs on a 250 x 250 size 
lattice at T = 250K and 6 — 0.1 ML. 

The general shape of the distributions seems to resem- 
ble previous results |35| for i* — 2,3 where the peaks 
rise more slowly on the left hand size (small s/s) due 
to the depletion of the mobile islands. This trend is 
qualitatively similar to the results for the case where 
small islands are unstable, as shown in Ref. |50|,[34|]. Note 
also that the peak height increases considerably as F de- 
creases. This may be due to coalescence which is found 
to become more pronounced as the deposition rate de- 
creases. Coalescence causes s to increase, pushing up the 
scaled island size distribution which includes the factor 
s 2 /9. 



A. Small-Island Mobility 

For both the EAM barriers and models I and II we ob- 
tain significant dimer mobility. The trimer mobility for 
the EAM barriers is rather small, and somewhat larger 
in models I and II. In general, we find that the energy 
barriers relevant for dimer and trimer mobility are in the 
same peak as the single atom hopping. The differences 
in hopping rates between monomers, dimers and trimers 
are due to very small differences in the energy barriers, 
which at T = 250K are significant. These differences in 
hopping rates decrease as the temperature increases. In 
model II, which has only one activation barrier for each 
peak, the activation barriers for monomer, dimer and 
trimer mobilities are all equal. Combinatorial summa- 
tion of paths then shows that for model II the diffusion 
coefficient of dimers (trimers) is equal to one half (one 
quarter) of the monomer diffusion coefficient, at all tem- 
peratures. Therefore, in model II, at any temperature 
in which atoms are mobile, dimers and trimers are also 
mobile. The mobility of small islands has a significant 
effect on the asymptotic scaling of the island density vs. 
deposition rate. In principle, this can be used to extract 
these diffusion properties at the atomic scale from STM 
and diffraction results at larger scales. However, our sim- 
ulations indicate that for experimentally relevant param- 
eters the system is likely to be away from the asymptotic 
regime, and therefore one should be careful in drawing 
conclusions from the empirically determined scaling ex- 
ponents. The relation between mobility of small islands 
and edge mobility was examined in Ref. J35|. It was 
shown that both types of mobilities are generated by the 
same hopping moves and are therefore related, namely, 
small island mobility implies edge mobility. 



B. Edge Mobility 

We have found that for Cu(001) atom mobility along 
straight island edges (for which the barriers belong to 
peak II) is much higher than the single atom mobility 
(with a barrier in peak III). This indicates, using an ar- 
gument from Ref. |g8j , that islands should form compact 
shapes for a broad range of temperatures. This conclu- 
sion is confirmed by our MC simulations and is in agree- 
ment with STM results ||. For the temperature studied 
here detachment of atoms from islands is negligible and 
therefore edge mobility is the dominant process which 
shapes the islands. 



IV. DISCUSSION 



C. Vacancy Mobility 

It is found that for Cu(001) the mobility of a single va- 
cancy is higher than the single adatom mobility. We ob- 
serve that the barrier for the diffusion of a single vacancy 
is sensitive to the environment beyond the 3x3 square 
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used to calculate energies in this paper. The energy bar- 
rier for diffusion of a single vacancy is E 1 ^ 1 = 0.12 eV 
if all sites around the 3x3 square are vacant. Adding 
three more atoms on the left it increases it to 0.28 eV. 
Embedding the hopping vacancy in an occupied 5x5 
square increases the barrier to 0.31 eV while the barrier 
for a single vacancy in the substrate is found to be 0.43 
eV (using a slab of 20 layers, 121 atoms/layer). We con- 
clude that, although the vacancy mobility is higher than 
the adatom mobility, at high coverage the difference may 
not be dramatic j37|. For example, at low temperature a 
difference of 0.05 eV in the barrier may significantly af- 
fect the hopping rate. A more general conclusion is that 
models I and II best describe diffusion and nucleation at 
relatively low coverage. At high coverage, far beyond the 
percolation threshold hopping can be considered primar- 
ily as motion of vacancies. When the adlayer is already 
dense, relaxation effects during the moves are more im- 
portant than at low coverage, and the configuration of 
occupied sites beyond the 3x3 square used here may 
affect the energy barriers. In spite of the prediction that 
a single vacancy is more mobile than an adatom, note 
that for Ag(001) there is evidence that vacancy clusters 
have lower diffusivity than islands [56| ], 

D. Comparison with Experiments 

The simulations presented here were carried out with 
parameters typical to homoepitaxial growth experiments 
on Cu(001) surfaces, namely, linear lattice size of 250 
sites (about 650 A), T = 250K and F between 4 x 10~ 4 
and lO^ 1 ML/s @|[5||. In Ref. ||. 7 = 0.33 was mea- 
sured for T = 223K, F between 2.5 x 10~ 4 and 10~ 3 
ML/s, and coverage 9 = 0.3. This value of 7 is in a very 
good agreement with our simulation results. The authors 
of Ref. jl5| claimed that this indicates that in this system 
and parameter range, i* = 1. Our results indicate that in 
this parameter range the system is away from the asymp- 
totic regime. Thus, the value 7 = 0.33 can be obtained 
eventhough i* — 2 or 3. A value of 7 = 0.55 was reported 
(for T = 223K, and F between 10~ 3 and 5 x 10~ 3 ) and 
7 = 0.58 (for T = 263K, and F between 2.5 x 10~ 4 and 
10~ 3 ) in the same coverage, 9 — 0.3 fli5| - These values 
are close to those reported in Ref. [jy], namely 7 = 0.54 
(for T = 220K, and F between 6.7 x 10" 4 and 3.3 x 10~ 3 ), 
and 7 = 0.46 (for T = 230K, and F between 4 x 10~ 4 and 
10~ 2 ). A possible interpretation of the higher values of 7 
observed in these experiments, is that processes such as 
coalescence become significant as the coverage increases. 

E. Diffusion on Other FCC (001) Metal Surfaces 

We expect our approach to apply also to other 
FCC(001) metal surfaces such as Ni(001) and Ag(001) in 
which diffusion occurs through hopping rather than by 



the exchange mechanism. EAM calculations for Ni(001) 
and Ag(001) are consistent with our models [^tJ. For sur- 
faces in which the exchange mechanism is favorable, as it 
is believed to be the case for Al @, Pt @, Pd and Au 
p9| modifications of the models are required. However, 
we believe that the approach presented here should still 
provide a simple and useful model for diffusion on these 
surfaces. 



V. SUMMARY 

We have studied adatom self diffusion and island 
growth on Cu(001) using MC simulations at the atomic 
scale. As input to the MC simulation we used a complete 
set of energy barriers obtained from the embedded atom 
method. To reduce the number of parameters and to ob- 
tain better understanding of the diffusion processes we 
have constructed two simple models. Model I, with three 
parameters, provides a good quantitative description of 
the full landscape of hopping energy barriers. Model II, 
with only two parameters, is a minimal model which still 
incorporates the essential features of adatom diffusion on 
the Cu(001) surface. 

We examined the diffusion properties of small islands 
and found that the mobility of dimers is comparable to 
the single adatom mobility while trimers are also some- 
what mobile. The mobility of adatoms along island edges 
was found to be much higher than the single adatom mo- 
bility. Since atoms detachment from islands is negligible, 
edge mobility is the dominant process that shapes is- 
lands into the compact forms. A further conclusion from 
the EAM calculations is that vacancy mobility, which is 
dominant at very high coverages, is higher than the single 
adatom mobility, which is dominant at low and interme- 
diate coverages. 

MC simulations of island growth show similar mor- 
phologies for the EAM barriers and models I and II. In 
all cases the islands generally form square or rectangu- 
lar shapes with a small number of kinks. Studies of the 
scaling of the island density N vs. deposition rate F 
show that N ~ F~< where 7 = 0.32 ± 0.01 for the EAM, 
0.33 ±0.01 for model I and 0.31 ±0.01 for model II barri- 
ers. These results are lower than the asymptotic value of 
7 > 0.4 anticipated for systems that exhibit dimer (and 
even trimer) mobility. It indicates that for the experi- 
mentally relevant conditions studied here the system has 
not reached the asymptotic regime. This result call for 
caution in using scaling results to draw conclusions about 
diffusion properties at the atomic scale. 

Since both models I and II have only few parameters, 
one needs a very small set of calculated energy barriers 
to determine these parameters for a given FCC(001) sub- 
strate, and a few more barriers to verify that the model 
applies to that substrate. This may open the way for a 
more effective use of first principle calculations of energy 
barriers on the surface as an input to the kinetics calcu- 
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lations. We believe that the procedure described in this 
paper would prove useful much beyond the FCC(OOl) 
metal substrates. We expect it to be applicable with 
some modifications to take into account exchange moves 
which are believed to be favorable in Al,Pt,Pd and Au. 
We believe that the approach presented here will provide 
useful models for adatom self diffusion on FCC(lll) sur- 
faces, for which reliable calculations of energy barriers 
are becoming available |6^j6l[], and also to a variety of 
heteroepitaxial systems. 
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Table Captions 

TABLE I. The hopping energy barriers obtained from the 
EAM calculations for all 128 possible configurations within 
a 3 x 3 square around the hopping atom. The barriers are 
given in eV. Each number in the Table is the barrier E B for 
the configuration in which the occupied sites are the union of 
the occupied sites in the picture on top of the given column 
(indexed by ni) and on the left hand side of the given row 
(indexed by U2). Consequently, the index n specifying the 
barrier is given by n — n\ + ni . 



TABLE II. Diffusion coefficients of small islands of sizes 
s = 1,2,3,5,7 atoms as measured from MC simulations of 
single islands. Results are shown for the EAM, model I and 
model II barriers at T — 250K. Islands of sizes s — 4, 6 and 
s > 8 were found to be immobile at this temperature. 
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Table II 



Cluster 




Diffusion Coefficients [hops/s] 




Size 


EAM 


model 1 


model 2 


1 


172 ±5 


118 ±4 


24.3 ± 0.8 


2 


237 ±8 


55 ±2 


13.2 ±0.4 


3 


5.6 ±0.2 


12.4 ±0.4 


5.9 ±0.2 


5 


1.12 ±0.04 


0.94 ±0.03 


1.00 ±0.03 


7 


0.28 ±0.02 


0.34 ±0.01 


0.42 ±0.02 
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Figure Captions 



Fig. 1: Classification of all possible local environments of a hopping atom including seven adjacent sites. Each 
site can be either occupied or unoccupied, giving rise to 2 7 = 128 local environments. Sites 1, 3 and 5 are nearest 
neighbors of the original site while sites 1, 2, 5 and 6 are adjacent to the bridge site that the atom has to pass. 

Fig. 2: The distribution of hopping energy barriers obtained from EAM calculations (a) and from model I (b) for 
the 128 possible local environments of Fig. 1. According to model I the lowest and highest peaks include 16 barriers 
each, while the middle peaks include 48 barriers each. The EAM peaks are broad and some overlap occurs. 

Fig. 3: Snapshots of surface morphologies obtained in MC simulations using the EAM (a), model I (b), and model 
II (c) barriers. The simulations were carried out on a 250 x 250 lattice at T = 250K and deposition rate F = 0.01 
ML/s. The snapshots presented are at coverage 9 = 0.1 ML. In all three cases the islands form generally compact 
square and rectangular shapes with a small number of kinks. 

Fig. 4: The island density vs. coverage is presented for the EAM (*), model I (o) and model II (+) barriers. In all 
cases the island density quickly increases at low coverage. It then saturates and remains nearly constant within the 
quasi-steady state aggregation dominated regime. The solid line is a guide to the eye. 

Fig. 5: Snapshots of the island morphology obtained from MC simulations using the EAM barriers for three 
deposition rates: (a) 0.1; (b) 0.01 and (c) 0.001 ML/s. The simulations were done on a 250 x 250 lattice at T = 250K 
and the snapshots were taken at 9 = 0.1 ML. It is observed that decreasing the deposition rate gives rise to fewer and 
larger islands. 

Fig. 6: The island density vs. coverage is presented for the EAM barriers, at deposition rates (from top to bottom) 
F = 0.1, 0.04, 0.01, 0.004, 0.001 and 0.0004 ML/s. The island density decreases as the deposition rate is lowered. Also, 
the plateau associated with the aggregation dominated regime broadens. The solid line is a guide to the eye. 

Fig. 7: The island density vs. deposition rate for the EAM (*), model I (o) and model II (+) barriers. These 
results were obtained at T — 250K and 9 = 0.1 ML for deposition rates between F = 0.1 and 0.0004 ML/s. The best 
fit through the six data points in each of the three models gives 7 = 0.32 ± 0.01 for the EAM, 0.33 ± 0.01 for model 
I and 0.31 ± 0.01 for model II barriers. These results indicate that models I and II provide a good description of the 
scaling properties. 

Fig. 8: Scaled island size distributions are shown for the EAM (a), model I (b), and model II (c) barriers, s is 
the average island size. The deposition rates examined in each case are F = 10 _1 (*), 10~ 2 (o) and 10~ 3 ML/s (+). 
These results are based on 20 runs on a 250 x 250 size lattice at T = 250K and 9 = 0.1 ML. 
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